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I We present a systematic survey of scarring and symmetry effects in the stadium billiard. The 

\l ■ localization of individual eigenfunctions in Husimi phase space is studied first, and it is demonstrated 

that on average there is more localization than can be accounted for on the basis of random-matrix 
theory, even after removal of bouncing-ball states and visible scars. A major point of the paper is that 
■ symmetry considerations, including parity and time-reversal symmetries, enter to influence the total 

' amount of localization. The properties of the local density of states spectrum are also investigated, 

as a function of phase space location. Aside from the bouncing-ball region of phase space, excess 
. localization of the spectrum is found on short periodic orbits and along certain symmetry-related 

lines; the origin of all these sources of localization is discussed quantitatively and comparison is made 
with analytical predictions. Scarring is observed to be present in all the energy ranges considered. In 
light of these results the excess localization in individual eigenstates is interpreted as being primarily 
due to symmetry effects; another source of excess localization, scarring by multiple unstable periodic 
^ ■ orbits, is smaller by a factor of 

(N 

(3 ■ I. INTRODUCTION 

o 
o 
o 



According to scar theory, the quantum eigenfunctions of a classically chaotic dynamical system do not always 
look locally like random superpositions of plane waves with fixed energy, as predicted by Berry instead, many 
eigenfunctions display a concentration of amplitude around short unstable periodic orbits greater than that expected 
on the basis of random matrix theory fluctuations. The first examples of scarring in the stadium billiard were presented 
by Heller Q , who also gave a semiclassical theory of scarring based on dynamics in the linearizable region around the 
periodic orbit (see also unpublished numerical work by McDonald Recent developments [Q have extended the 
^ . theory of scars to the non-linearizable regime, to include the effects of homoclinic recurrences at long times. Scarring is 
then seen to be a weak localization phenomenon coming from the short-time correlations associated with an unstable 
periodic orbit; in the energy domain, this means that a wavepacket centered on the orbit may have non-random 
overlaps with the eigenstates of the system. Quantitative measures of the strength of scarring have been developed 
and tested numerically Q,^. 

In this paper we ask whether the localization caused by scarring on not-too- unstable (i.e., short) periodic orbits, 
and by atypical regions such as the "bouncing ball" modes is adequate to predict measures of localization in 

eigenstates of chaotic systems. Quantitative numerical confirmation of scar theory jj,^ was limited to discrete-time 
maps, and it would be desirable to study scarring quantitatively in the context of more realistic and experimentally 
realizable systems, such as the stadium (Bunimovich) billiard. In what follows we present a systematic study of 
eigenstate localization in the stadium billiard (see Ref . ||] for some other recent analyses of this system) . In Section ||, 



X 



we first present the numerical method used to find the eigenstates. Then, in Section [II, we study the localization 



properties of individual eigenstates, followed by an investigation in Section [V of the localization properties of the local 
density of states. We find evidence of scarring as a ubiquitous phenomenon, in all the energy ranges considered. We 
examine quantitative predictions of scarring strengths based on the classical structure around the unstable periodic 
orbits Q. However, scar theory is not sufficient to explain all the observed localization. Symmetry effects, including 
the imprints of both time-reversal and parity symmetries, are ultimately found to dominate the excess wavefunction 
localization in the stadium. 
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II. METHOD 



We study eigenstates of the time-independent Schrodinger equation in an infinite 2-diniensional potential well in 
the shape of a stadium, taken here to consist of a square of side 2 with semicircular endcaps of radius 1. We use the 
plane wave method pO| to find states with even parity with respect to reflection about each of the two symmetry 
axes of the stadium. It is assumed that eigenstates at energy E can be approximated as a superposition of plane 
waves at that energy; i.e., with |kp = E (we use h = 1 and to = 1/2 here and throughout). Therefore we use as 
a basis plane waves with even-even parity, namely cos{kx cosOj) cos{ky sinOj), with N values of 9j chosen uniformly 
between and tt/2. Then the coefficients of these N plane waves are determined by minimizing the squared value of 
the wavefunction at M > N equally spaced points on the boundary. The wavefunction is set to an arbitrary non-zero 
value at a point in the interior in order not to have a singular system of equations, and M is kept slightly larger than 
N so that the system of equations is overdetermined. In practice M must be chosen so that there are several points 
per wavelength. 

Since the eigenvalues are not known a priori, we search for them using the 'tension' method: at each energy E we 
solve the linear system for the coefficients of the plane waves by singular- value decomposition; then, we compute the 
integral of along the boundary. If E is an eigenvalue this integral should vanish, to within numerical precision. 
In practice one finds sharp minima in the tension as a function of k that are identified with the eigenvalues. Since it 
is not known in advance how many plane waves are needed to give an accurate representation of the wavefunction at 
a given energy, we repeat the search for minima with increasing numbers of plane waves N and boundary points M 
until the results of subsequent iterations agree; we find, for instance, that for k around 100, N — 300 and M = 310 are 
necessary. This corresponds to eight points per wavelength along the part of the boundary lying in the first quadrant. 
The number of plane waves needed scales linearly with k. 

The computed eigenvalues and eigenfunctions satisfy three diagnostic tests: (1) the total density of states agrees 
well with the Weyl area rule; the deviation is never more than about 5 states out of 5000, or 0.1% (the bulk of 
this deviation is due to a small periodic modulation in the level staircase function, which has the right period to be 
attributable to the bouncing-ball states Q ; after taking out this modulation we find that the total number of possible 
missing or extra states in the region 100 < k < 200 is not more than one or two); (2) the histogram of level spacings 
is consistent with the Wigner surmise for random-matrix statistics; (3) the overlaps between eigenstates are found to 
be less than 1%. A problem with our procedure is that it produces many shallow tension minima, by which we mean 
that the contrast between the tension at a local minimum and the background tension is less than a factor of 100. 
The frequency of such shallow minima grows with increasing energy. Moreover, these poor minima are not improved 
by increasing the parameters M and N because doing so makes the matrix more ill-conditioned. This suggests that 
there is an intrinsic shortcoming in the ability of our basis to represent eigenstates; to have a perfectly accurate 
representation one would presumably need to include evanescent waves of total energy E as well. Although it has 
been shown by Berry ||Tl| ] that an evanescent wave may be represented by a singular superposition of plane waves, 
this representation is numerically ill behaved, giving poor convergence properties. 

In spite of this difficulty, essentially all of the questionable states must correspond to true eigenstates, in order for 
the counting of states to come out right. Also, many of the shallow minima we assume correspond to eigenstates 
occur in close conjunction with other minima, and yet the numerically computed histogram of level spacings decreases 
linearly toward zero at small spacings, in good agreement with the theory. 

Another difficulty with the present method is that the computation time per eigenstate scales as k^, making the 
collection of extensive eigenstate statistics above k — 300 almost prohibitively expensive. The plane-wave method 
described here has been improved upon by Li and Hu |l^ , who expand the derivative of the tension analytically in 
k in order to predict the position of the next minimum. Li and Hu's approach improves the speed of the plane-wave 
method by a factor of about five, but does not change the way the computation time scales with k. Thus, only slightly 
higher values of k could be reached for a given computation time, and we have not needed to employ this approach 
here. Another approach to the stadium billiard has been outlined by Saraceno's group ]l3[ |. Their method, suitable 
for very high energies, obtains all the eigenvalues in a narrow range of k simultaneously by solving a generalized 
eigenvalue problem, rather than searching for them step by step. The energies to be studied in this paper, however, 
are not so high as to necessitate the use of such an algorithm. 

Lists of eigenstates of the even-even symmetry class in the three energy ranges of fc = 20 to 30, k — 100 to 150 
and k = 200 to 225 were generated according to the method outlined above. Some examples are given in Fig. |l| of 
eigenfunctions ^'„(x,?/) in coordinate space. 
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A. Surface wavefunctions 



Classically, one typically uses the boundary of the billiard as a Poincare surface of section. The variables parame- 
terizing the surface are q, the arclength along the boundary, and its conjugate momentum p, which is the component 
of momentum parallel to the wall (both positive in the clockwise sense), q is taken modulo the length L = 4 + 2tt 
of the perimeter, while p is limited by energy considerations to —k < p < k. Then the classical billiard dynamics is 
reduced to a nonlinear one-bounce map, = Mci{q,p)'^ ■ 

A natural way of reducing the quantum problem from two to one dimension is to characterize the Dirichlet eigen- 
functions y), which are defined on the interior of the stadium, by their normal derivatives (pniq) on the boundary: 

=<?!)„ (g) cxn-V^'„(x,?/). (1) 

We will study the properties of these surface wavefunctions. The wavefunctions are normalized according to the 
convention 



dxdy\^„{x,y)\\ (2) 
which is equivalent to the following condition on the surface wavefunctions: 

^Jdq{h{q)-r{q))\Mq)\^^^, (3) 

where n(g) is the unit normal at the position q along the boundary, and r{q) is the displacement vector from the 
center of the stadium to the position q. 

The Husimi representation of a surface wavefunction \n) is given by its projection | {qo,po\n) p onto test state 
Gaussians of the form 



{q\qo:Po) = ^TT^TTiexp 



(g - g o) 

2^2 



2 



(4) 



The Gaussian is centered at the point {qo,Po) in the boundary phase space. The parameter a controls the aspect ratio 
of the Gaussian: its width in position is cr/\/2 while its width in momentum space is l/{a^/2). To maintain a given 
aspect ratio of the Gaussian in phase-space, a must scale as A:"^/^. Here, we choose a = [(2 + 7r)/fc]^^^, which makes 
the aspect ratio unity in units where the full phase space on the billiard boundary is taken to be a square. Note that 
the wavefunctions are real in coordinate space but complex in phase space. Some example Husimi plots are given in 
Fig. I 

III. IPR STATISTICS FOR INDIVIDUAL EIGENSTATES 

We start by investigating the localization properties of individual eigenstates. As a measure of the degree of 
localization of the eigenstates in phase space, we use the mean squared value of the intensity, also known as the 
inverse participation ratio (IPR) 

IPR„ = ' ' , , (5) 

where qj and pj range over the entire phase space as j varies on a sufficiently fine grid, which we take to be 200 x 200. 
The range of momenta pj covered is from —k to +k for a wavefunction with energy fc^; outside this classically allowed 
region there is almost no wavefunction amplitude. The IPR measures how much fluctuation across phase space there 
is in the eigenfunction. If \{qj,pj\n)\^ were completely uniform over phase space, the IPR would equal 1. On the 
other hand, if all the intensity were concentrated entirely at one point j and was zero elsewhere, the IPR would reach 
its maximal value J. [Of course, the uncertainty principle prevents the IPR from ever becoming greater than the size 
of phase space in units of a Planck cell, i.e., 0(fc_L/27r).] Random-matrix theory (RMT) predicts a Porter-Thomas 
distribution of wavefunction intensities, which yields an IPR of 2 (for complex wavefunctions), already larger than 
the naive classical expectation of 1. 
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A subtle point is that the wavefunction intensity |n) | , even when averaged over wavefunctions \n), is not 

uniform over the boundary phase space, being instead a non-trivial function of pj (see Eq. (|^) below). Therefore, 
even in the absence of quantum fluctuations, if the intensities \(qj ,pj\n)\^ for a wavefunction \n) were given simply 
by the mean value ^/l~^~(pJky as in Eq. |8[ we would obtain an IPR of 

(l-(p/fc)2)/(Vl-b/fc)^)'«l-08. (6) 

The averages (• • •) here are of course taken over the classically allowed range of p: ~k < p < k. If the intensity 
fluctuations around the smooth value -y/l — {p/k)"^ were instead given by a Porter-Thomas distribution, we would 
obtain an IPR of 1.08 x 2 = 2.16 (since the smooth behavior and the Porter-Thomas fluctuations are independent of 
one another, the IPR contributions simply multiply). 

In Table | the IPRs of the first 36 eigenfunctions in each energy range are tabulated. The average IPRs for the first 
100 states in each energy regime are summarized in Tabic ||. The IPRs cluster around 3, with some much larger. 
The mean IPR is 3.67, 3.60 and 3.08 in the low, medium and high energy ranges respectively, but with large standard 
deviations. Visual examination of the Husimi plots for states with the largest IPRs reveals that they correspond to 
bouncing-ball states. Since this phenomenon is well understood (and is associated with marginally stable periodic 
orbits), we remove this effect by removing the bouncing-ball states from the average (see column (e) of Table ||). 

The remaining states avoid the bouncing-ball region, causing their IPR to be increased by an amount that can be 
estimated. At fc « 20, the bouncing-ball region occupies about 13% of phase space. Assuming the non-bouncing-ball 
states to vanish identically in that region, the random-matrix theory prediction for the average IPR without the 
bouncing-ball states becomes 2.16 x 1.13 = 2.44. The size of the bouncing-ball region in phase space is energy- 
dependent (and must of course go to zero in the fc — > c» limit, in order to satisfy the Shnirelman quantum ergodicity 
condition). As found by Backer et al. the number of bouncing-ball states scales as NtbiE) = 0.2£'^/"' and thus the 
density of bouncing-ball states, which is proportional to the area occupied by them in phase space, scales as E^^^'^ 
or This scaling law is in agreement with our numerical results. The adjusted RMT predictions are shown in 

column (d) of Table 0. 

As seen in Table |j the average IPRs with bouncing-ball states removed are still well above the random-matrix 
prediction in each energy range. The IPRs were obtained by averaging over 100 states; the statistical uncertainty in 
the means is thus about l/VlOO « 0.1 of the standard deviation, or roughly 0.1, which is comparable to the differences 
in average IPR between different energy ranges. We see that the IPR does not increase markedly with wavevector k, 
so the stadium is not an example of weak quantum ergodicity as studied by Kaplan and Heller |]15[ (this makes sense, 
since the stadium is strongly chaotic with rapid mixing throughout its phase space with the exception of the bouncing 
ball region). However, neither is there a marked trend towards the RMT predictions as k increases by a factor of ten. 

Thus, we have uncovered systematic evidence of localization of eigenfunctions in the stadium billiard beyond what 
would be predicted on the basis of random-matrix theory. What is responsible for it? The next highest IPRs, after 
the bouncing-ball states, are for eigenfunctions visibly scarred along the periodic orbit that bounces horizontally 
between the centers of the endcaps (these "horizontal-bounce states" are marked with a dagger in Table ||). Thus, one 
might hypothesize that scarring along this and other periodic orbits is partly responsible for the excess localization. 
In the next Section, however, we shall introduce a method that allows the excess localization to be identified with 
specific classical structures in phase space, and techniques to predict that localization theoretically. We shall find 
that quantum symmetry effects cause most of the excess localization, while a secondary effect consists of combined 
scarring contributions from all of the short unstable periodic orbits. 

IV. SPECTRAL LOCALIZATION 
A. Introduction 

The IPRs of individual eigenstates establish that there is localization compared with the predictions of random 
matrix theory. However, it does not permit the localization to be definitively identified with structures in classical 
phase space because the IPR is an aggregate over all of phase space. But it is also possible to ask a complementary 
question, namely, how ergodic are the eigenstates at a given point in phase space as one varies the energy? To answer 
this question, we turn to the local IPR — a tool that provides a picture of the regions of phase space that are prone to 
intensity enhancement or depletion over an ensemble of eigenstates. 
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TABLE I. Inverse participation ratios of the Husimi plots of the first 36 eigenstates in the low {k = 20), medium (fc — 100) 
and high (fc = 200) energy ranges. Bouncing-ball states are labeled with an asterisk, and states visibly scarred on the 
horizontal-bounce orbit are labeled with a dagger. The random-matrix theory prediction, with bouncing-ball correction, is 
2.44, 2.29 and 2.25, respectively, in the three energy ranges. 
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TABLE II. Averages of the inverse participation ratios of the first 100 states in energy ranges starting at fc = 20, 100 and 
200. Columns: (a) the k value at the beginning of the range; (b) the actual IPR; (c) the estimated fraction of phase space 
occupied by the bouncing-ball states in this range of k; (d) the predicted IPR when the bouncing-ball region is excluded; (e) 
the computed average IPR after the bouncing-ball states (e.g., those marked with an asterisk in Table |) are omitted. Standard 
deviations are shown in parentheses. The large remaining discrepancy between (d) and (e) is mostly due to symmetry effects, 
as explained in the text. 
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1. Definition 



The local IPR ('LIPR') is the mean squared eigenstate intensity at point (g,p), averaged over an ensemble of 
eigenstates fig]: 



LIPR(g,p) 



1 ^ I 



n=l 



1^1 2\ 
n=l ^ 



2 ' 



(7) 



where TV is the number of eigenstates being summed over. Typically the sum is taken over eigenstates in a small 
energy range around some central value of k. The Gaussian wavepackets \q,p) are adjusted to maintain constant 
aspect ratio in phase space as k changes (see discussion following Eq. (Q)); this keeps the classical dynamics fixed as 
the energy increases. Several LIPR pictures are shown in Fig. ||. 



2. Intuitive description 



The LIPR at the point {q^p) is a statistical property of the set of 'random' variables {{q,p\n) : n = 1 . . . iV}. The 
second moment of this quantity [which appears in the denominator of Eq. ^] is a smooth function of q and p, as 
will be shown in Sec. IV C; it is present as a normalization factor. The lowest moment that is nontrivial is the fourth 
moment, which appears in the numerator of Eq. (|^). The LIPR, then, measures the non-uniformity of the local density 
of states at a (fixed) position in phase space. A large LIPR indicates that a relatively small fraction of eigenstates 
have large overlaps with the test Gaussian, whereas a small LIPR indicates that the overlaps are distributed in a more 
'egalitarian' way. 

We expect the LIPR to be a sensitive indicator of the presence of scarring, because a wavepacket centered on a 
periodic orbit has a local density of states that oscillates with energy Eigenstates in the peak of the oscillation 
will, on average, have enhanced overlaps with the test state, while eigenstates in the trough will have suppressed 
overlaps. This nonuniformity will cause the LIPR to have peaks near periodic orbits. 

There is a useful interpretation of the LIPR — it is proportional to the long time average of the probability that 
the wavepacket evolved in time, has returned to its original location. A derivation and discussion of this 

interpretation are presented in Sec. IVG below. 

Since the stadium is strongly chaotic, classical trajectories explore all of phase space (except for a set of measure 
zero) . Therefore the classical return probability for long times is uniform, and the resulting naive classical prediction 
is that LIPR — 1 . This is indeed too naive. 

Instead, the 'null hypothesis' for a chaotic system is that its eigenfunctions are random superpositions of plane 
waves Under that maximally random assumption, projections of wavefunctions on test states should follow a 
distribution with one degree of freedom (if the eigenfunctions and test functions are real) or a distribution with 
two degrees of freedom (if they are complex). The corresponding predicted LIPRs are 3 or 2, respectively — already 
different from the naive classical prediction. 



B. Introduction to pictures 

In Fig. ^ we present the LIPR computed for even-even eigenstates in the wavenumber ranges fc = 50 to 60, 100 to 
150 and 200 to 225. They display interesting and beautiful localization effects. 

The horizontal axis of each picture is q, the position along the boundary, which runs from to L. q — Q corresponds 
to the midpoint of one of the straight walls. The vertical axis is the component of momentum parallel to the boundary 
p/k, which runs from —1 to 1. The eightfold symmetry of these plots results from the two spatial reflection operations 
plus time-reversal symmetry. The LIPRs all show a very large amplitude in the bouncing-ball region, explainable by 
the bouncing-ball and near-bouncing-ball states present in the energy ranges considered. But in addition, there is 
localization in many other regions of phase space. There are very prominent straight lines of enhanced LIPR running 
along the lines of symmetry, and filamentary streaks running across other parts of phase space. There are also isolated 
spots at which the LIPR is unusually large. 

Note that, although the details vary somewhat, the main localization features appear in the same positions in each 
energy range. This lack of ?i-dependence suggests that the streaks are associated with classical phenomena. Our 



6 



main goal is to give a quantitative explanation for all of these effects. We would also like to understand the relative 
contributions of each semiclassical effect to the total average IPR enhancement found in Sec. Ill above. 

If one varies the aspect ratio of the test state to make it more position-like or momentum-like, a series of bright 
spots can also be resolved at various points along the streaks, on a scale as small as the smaller axis of the Gaussian 
test state (see Fig. ^). 



C. Mean wavefunction intensity 

The mean wavefunction intensity appears in the denominator of Eq. (|^) in the role of a normalization factor. It can 
be determined from classical phase space arguments similar to those used in the derivation of Weyl's law. It depends 
on p for two reasons: (1) we are using normal derivatives, which introduces a contribution of p± — k^l — (p/k)'^ 
multiplying the coefficients of the plane waves (since the coefficients are squared, the resulting factor is fc^(l — (p/fc)^)), 
and (2) there is a geometrical factor coming from the projection from the circle k'^ + p'^ = k^ onto the boundary of 
the stadium. Because the plane waves are evenly distributed around the circle, there are more near p — k than near 

p = 0; the resulting geometrical factor is [l — (p/fc)^] . Together, (1) and (2) give 



[q,p\n) 



|2 



Vl-ip/k?. (8) 



This smooth factor is divided out of the LIPR in Eq. (^, giving the LIPR a fiat background. 

We checked the mean wavefunction intensity for the computed even-even states in the interval from k = 100 to 110. 
It agrees well with Eq. (^), except at p — around the centers of the straight segments and the centers of the endcaps, 
where sharp peaks are seen. The peaks occur because only even-even states have been included in the analysis; since 
the odd states must vanish at the symmetry points, the even-even states compensate by having double the average 
intensity there. Apart from this, the agreement of the computed mean intensity with Eq. (^) indicates that we have 
included enough eigenfunctions in the energy ranges chosen for Fig. ||. 



D. Symmetry lines 



What is the explanation for the streaks in the LIPR plots of Fig. g? The most prominent streaks are shown 
schematically in Fig. ^. The streak labeled (1) in Fig. ^ corresponds to trajectories emerging from the center of the 
endcap. The streak labeled (2) corresponds to a family of orbits emerging from the center of the straight segment. 
Streak (3) , which starts at small momentum at the edge of the bouncing-ball region and extends to large momentum 
at the center of the endcap, corresponds to orbits that have a vertical segment in the endcap region. The box and the 
bowtie are the two most prominent periodic orbits in this family. The streak labeled (4) starts at small momentum at 
the center of the straight segment, curves up to larger momentum, and finally comes back down to zero momentum at 
the center of the endcap. It corresponds to orbits that pass through the center of the stadium. The most prominent 
periodic orbit in this family is the Z orbit. Finally, there is the horizontal line (5) going across the plot at zero 
momentum, corresponding to normal incidence on the billiard boundary. All of these families can be understood by 
considering time reversal and parity symmetry effects in combination with dynamics, as we will now see. 



1. Simple explanation 

The intuitive reason for the symmetry lines is that the system, having time-reversal symmetry, has real eigenfunc- 
tions. However, the Gaussian test states [Eq. (^] are intrinsically complex. Therefore the overlap {q,p\n) has both 
real and imaginary parts which, for most choices of q and p, are not trivially related to one another. Therefore the 
typical computed LIPR for the system is 2, characteristic of complex eigenfunctions. But in a sense the 'correct' value 
LIPR = 3 (characteristic of real eigenfunctions) is obscured by a non-ideal choice of test states whose properties do 
not match those of the eigenstates. 

However, due to the symmetry of the system, there are certain values of q and p for which the real and imaginary 
parts of {q,p\n) are not random and independent. For example, when p — (normal incidence to the wall), the 
test states and their projections {q,p = 0\n) on the eigenstates become pure real and therefore near p = the LIPR 
increases to 3. [This explains the symmetry line labeled (5) in Fig. |^.] Similarly, the eigenfunctions are symmetric 
about the centers of the straight walls, and therefore {q = 0,p\n) is pure real and again the LIPR increases to 3. 
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(Similarly for the centers of the endcaps.) This simple argument explains the straight symmetry lines (1), (2) and (5) 
that appear in the LIPR pictures. A rigorous derivation of their heights and widths will be given in the next Section. 

It is interesting to note that the same LIPR enhancements would be present if we had used eigenstates of other 
symmetry classes, or put eigenstates of all symmetry classes together. For example, for odd-odd eigenstates, the 
eigenf unctions are antisymmetric about the centers of the straight walls, and therefore {q = 0,p\n) is pure imaginary, 
again giving a distribution with one degree of freedom, and the enhancement LIPR = 3. 



2. Derivation 



In this Section we show how to compute the enhancement of the LIPR that appears near the symmetry lines of the 
system. We do this in order to verify the conclusions of the preceding simple arguments, and also to deduce important 
details such as the profiles and widths of the symmetry lines. 

The derivation relies on the fact that the eigcnfunctions are chosen to be be real (in a coordinate basis), which can 
be done as a consequence of time-reversal symmetry. We denote time-reversal with T, an antiunitary operator. On 
the surface of section, T\q,p) = \q, —p). 

Furthermore, the eigcnfunctions are symmetric or antisymmetric with respect to the reflections x — > —x and 
y —y, which we will denote with the unitary operators Rx and Ry respectively. On the surface of section, this 
corresponds to symmetry with respect to reflection about two values {qx, Qy} corresponding to an intersection of the 
billiard wall with the y and x axes, respectively. Then for example, Rx \qx + A(7,p) = \qx — Ag, —p) (where q is always 
taken modulo the billiard perimeter). For simplicity in the following derivation we will consider the simpler case of a 
single reflection symmetry, R, about q — 0. 

It is not correct to model the eigcnfunctions as uncorrelated Gaussian random variables, because their reflection 
symmetry correlates their values at different points. However, we can generate random wavefunctions |n, ±) with 
positive ('-!-') or negative ('— ') symmetry by taking completely random real wavefunctions |n) (with no symmetry) 
and projecting them onto the correct symmetry subspace: 

\n,±)^^{\n)±R\n)). (9) 

To substitute into Eq. (^, we will need to compute {q,p\n,zL). This quantity is complex, so we decompose it into 
two real random variables ^± and i'± as 



1 

72 



— [{<Pq,p\'n) ± (09,pl^l")] = ^J■±+ iv±. (10) 



If /i± and i'± had the same variances and were uncorrelated, the LIPR would uniformly equal 2, like that of any 
Gaussian random process with two degrees of freedom. This null result would hold even if the variance depended on 
p and g, because the square of the variance appears both in the numerator and in the denominator of the definition 
of the LIPR. 

However, in reality the symmetries cause the variances of the real and imaginary parts ^± and v± to depend 
differently on phase space location. The effective number of degrees of freedom varies from two (when (^j.) — {v'±)) 
to one (when (/ij.) ^ 0, (i^^.) — 0), and correspondingly the LIPR varies from 2 to 3. 

In terms of the three quantities (/Xj_), (i^^), and (/i±i^±), 

(|m± + *^±|') = (4) + (4), (11) 

and ^|/i± -l- iv±\'^^ can be expanded and then contracted pairwise, giving 

(|a*± + *^±l') = 2 ((4) + + ((4) - + 4(a.±^±)^ (12) 

The problem is reduced to the computation of the required variances and correlations in Eq. (p^). 
Using {q,p\n)* — {q, --p\n) and R \q,p) = \—q, — p), we obtain 

M± = [{q,P\n) + {q, -p\n) ± {-q, -p\n) ± {-q,p\n)] (13) 

= [{q,p\n) - {q, -p\n) ± (-(?, -p\n) =F {-q,p\n)] . (14) 
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Because \n) represents not a wavefunction but rather the projection of a wavefunctio n onto t he surface of section, the 
closure relationship does not hold: ^„ |"-)("'| 1- But as explained above in Sec. IV A 1 , it is still approximately 
true that 



X! I")("k'-P> ~ lip) k'-p) ' 



(15) 



where f{p) = ,f{—p) ^ \/l — {p/ky. It follows that 

n 

X [{n\q,p) + {n\q, -p) ± {n\-q, -p) ± {n\-q,p) 
{q,p\q,p) +^e{q,p\q,-p) ±Re{q,p\~q,p) ± {q,p\-q, -p) 



fip) 



2N 



(16a) 



the last line follows from {q, —p\—q,p) — {q,p\—q,—p), which is real. The other correlations can be worked out 
similarly; the results are 



(t^i) = {q,p\q,p) -'l^e{q,p\q,~p) T^e{q,p\-q,p) ± {q,p\-q,-p) 

f {_P) ' 

(m±J^±) = Im((7,p|g,-p) ±Ini((7,p|-q,p) 



(16b) 
(16c) 



3. Even and odd subsets vs. whole set 

At this point we must distinguish between two types of LIPR. First, one can evaluate the LIPR by averaging (in 
the numerator and the denominator) over only the even states (as done in this paper) or only the odd states. In that 
case one obtains 

LIPR(g,p; {|n, ±)}) = 2 + K^Ml, -p) ± {lM-1,p)\'^ (^^^^ 

i{q,p\q,p) ± {q,p\-q,-p)) 
^ 2 + l(g^Pl^lg'P) ^ {Q,P\RT\q,p)f ^^^^^ 
{{q,p\q,p) ± {q,p\R\q,p)f 

On the other hand, one can evaluate the LIPR by averaging over both the even and the odd states. In that case, 

TTPT^^ /I _l\\mji o , \{q,p\q,-p)f + \{q,p\-q,p)f x 

LIPR(9,p;{|n,+)}U{|n,-)}) =2 + — ■ ■ (18a) 

{q,p\q,p) + {q,p\-q,-p) 



\{q,p\T\q,p)\^ + \{q,p\RT\q,p)\^ 

{q,p\q,pf + {q,p\R\q,pf 



The required matrix elements can be worked out from Eq. (^: 

{q-q'f o\p-p'f 



\{q,p\q' ,p') \ = exp 



4ct2 



(18b) 



(19) 



The three LIPRs — Eq. (|17| ) for even or odd symmetry and Eq. ( |18D for both symmetry classes put together — all can 
be shown to give results between 2 and 3. In fact, outside the region near q = p = Q the three versions are identical; 
they all predict bright symmetry lines near q — Q (or near any q associated with a parity symmetry) of the form 

LIPR(q, \p\ > l/a) = 2 + 6-2-?'/^' (20) 

and bright symmetry lines near p = of the form 

LIPR(|q| >CT,p) = 2 + e-2p'-'. (21) 

We note that the point p = g = 0, near which the three LIPR definitions given above do not agree, will always be 
on a short periodic orbit (e.g., the horizontal bounce orbit of the stadium billiard), and so the be havior there in any 
case cannot be determined without taking dynamical scar effects into account (see Section IVE2). 
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4- Curved symmetry lines 



The curved streaks (3) and (4) also correspond to classical structures related to the system's symmetry. Consider 
first the streak labeled (3) in Fig. ^. It corresponds to trajectories with a vertical segment in the endcap region. In 
the desymmetrized quarter-stadium, these trajectories are incident normally on the y = horizontal boundary, so a 
Gaussian test state with zero tangential momentum placed at the point of normal incidence {q\p = 0) would give 
rise to an enhanced LIPR of 3, as for the zero-momentum line (5) above. Of course, this point of normal incidence 
exists only on the boundary of the gworter-billiard. At other points along the same orbit, including points that also 
live on the boundary of the full stadium, a Gaussian with momentum aligned along the trajectory will be close to, 
but not exactly equal to, a time iterate of this zero-momentum Gaussian living on the y = boundary. The exact 
time iterate of this purely real Gaussian will in general be a Gaussian centered on the other periodic point {q,p) but 
with a (possibly complex) width a somewhat different from that used in our test state. Thus, our test Gaussian at 
((7,p) will have significant overlap with an iterate of a purely real state at {q' ,p = 0) and result in an enhanced LIPR 
somewhere between 2 and 3. Very little wavepacket rotation or stretching occurs during the short vertical trip between 
the y = line and the endcap; it is for this reason that the streak (3) is so strong. Further iterations of the dynamics 
will result in additional streaks of enhanced IPR; these however will be much weaker due to the additional stretching 
which makes a circular wavepacket centered at these distant points less closely related to the evolved version of the 
{q' ,p = 0) real wavepacket. All streaks in the LIPR plots which have not been explicitly identified in Fig. ^ may be 
explained as further iterates of symmetry lines we have discussed explicitly. The rather strong streak (4), associated 
with trajectories going through the center of the stadium, has a completely analogous explanation. 



5. Quantitative contribution of the symmetry lines to the IPR 



We can now estimate the contribution of the symmetry lines to the predicted average IPR. There are twelve strong 
roughly vertical streaks (including the curved streaks (3) and (4)) in the region p 7^ and one horizontal streak at 
p = 0. We assume temporarily the curved symmetry lines to have a central height near 3 and width of order a just as 
for the a; = symmetry line. Integrating the area under these streaks assuming Eqs. ( ^ ) and (|l]) gives a predicted 
IPR from symmetry effects alone of 2.51 (in the energy range from k = 100 to 150). The enhancement of the LIPR 
in the whispering-gallery region and in bright spots due to periodic orbits gives a contribution only of about 0.03, 
negligible in comparison to the symmetry effects. Since the bouncing-ball region, which covers about 6% of phase 
space, is excluded in the above analysis, we should increase the prediction by 6% to 2.66 for the purpose of comparison 
with the numerical value of 2.70 ± 0.07. In the energy range from k = 200 to 225 the corresponding predicted IPR is 
2.40, or 2.51 corrected for the bouncing-ball region, whereas the numerical value is 2.82 ± 0.08 (the discrepancy in the 
higher-energy region remains unexplained). Thus, it seems that inclusion of symmetry and bouncing ball effects is 
sufficient to yield a quantitative understanding of much of the excess average IPR in Table ^ Interestingly, scarring 
plays little role in the excess average IPR, but a very important role in understanding the phase-space structure. 



6. LIPR using real wavepackets 



In light of our claim in Sec. |IVD l| that complex Gaussian test states living on the boundary of the billiard are not 
ideal for probing eigenstates in a system with time reversal symmetry, it is natural to ask what test states would be 
more appropriate. A natural approach would be to symmetrize the test wavepackets with respect to time reversal 
symmetry, by taking their real or imaginary parts. Equivalently, since the eigenstates are chosen to be pure real, 
one may simply take the real or imaginary parts of the overlaps {q^p\n). We note that one might also consider 
symmetrizing the test Gaussians with respect to parity symmetry, which is also present in our system. This, however, 
would have no effect on any observed quantities, since all of the eigenstates already respect parity symmetry. We then 
define 



LIPR,eal(g,p) 



1 ^ I 



p\n) 



1^1 2\ 
— ^|Re(g,p|n) j 



(22) 
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This quantity is plotted in Fig. |5[ and should be compared with the corresponding pictures in Fig. || for the LIPR of 
the original complex wavepackets. Because the new test states are real, the background LIPR rises from 2 to 3, while 
the LIPR along the straight symmetry lines remains unchanged (at 3). 

Naively, one might expect the desymmetrization process described here to completely eliminate all symmetry effects, 
leading to flat LIPRrcai('?,p) = 3 behavior, with the exception of dynamics-related enhancement in the bouncing-ball 
region and near isolated unstable periodic orbits. Instead, one finds that although the straight symmetry lines (p = 
and the parity symmetry lines) have indeed disappeared in Fig. ^, LIPRjoai along the curved symmetry lines (e.g. the 
streaks identified as (3) and (4) in Fig. U) is enhanced from 3 to ~ 4.5, making these curved streaks just as visible 



as in the original plot of Fig. g. This is easy to understand based on our analysis of these streaks in Sec. IV D 4 
We saw there that wavepackets located on the curved streaks are time-evolved versions of real wavepackets. So let 
r„ be the overlap of eigenstate \n) with a real wavepacket (i.e. either one having p = or one that lives on a parity 
symmetry line). r„ of course is distributed as a real Gaussian variable. Now the overlap of the same eigenstate \n) 
with a time-evolved copy of the real wavepacket will be given by = r„expi0„, where in the semiclassical regime 
6n can be taken to be a random phase. The LIPR for r'^^ is of course equal to that of r„ (i.e., 3), as the phase has 
no effect on the intensity = knP- Ifj however, we take the real part of the wavepacket living on the curved 

symmetry line, we must look at the quantity 

LIPRreal = /^r'^^'""^ = I (23) 



2 



because Re r'^ — r„ cos0„, where r„ is Gaussian-distributed and On is a random phase. This is in contrast with the 
LIPR of 3 obtained on the straight symmetry lines as well as in the background (the LIPR on the straight symmetry 
lines would also be enhanced to 9/2 if we had let the eigenstates have random phases instead of adopting the convention 
where they are all real). Thus, there appears to be no natural and simple way to eliminate all symmetry effects on 
the LIPR (and thus, also on the IPR) in this time-reversal invariant system. 

E. Dynamics 

In this Section we will consider the LIPR for a system with no time-reversal symmetry, where the wavefunctions 
are complex instead of pure real, and the arguments of Section IV D| concerning symmetry lines do not apply. In such 



a system the LIPR plot would be far more uniform, and most of the remaining structure would consist of isolated 
regions of enhanced LIPR, against a background of LIPR — 2 (corresponding to the statistics of the square of a 
complex Gaussian variable). The streaks would be all gone, including the vertical parity-line streaks, which also 
depend on time-reversal symmetry for their existence. Thus in this Section we study the effect of dynamics alone on 
the LIPR. First we give a general theory for computing the LIPR classically anywhere in phase space, and then we 
show an easier way to compute the LIPR in the neighborhood of a periodic orbit. In the process, we will also see how 
to treat the case of a p eriodi c orbit that happens to lie on a symmetry line in a time-reversal invariant system. In the 
following Section (Sec. [VF| ), we will describe an alternative way to treat parity symmetry and dynamical effects in 



a unified framework, to better understand scar effects in a system with symmetry. 

1. Full dynamics treatment 

In this Section we show how to compute the enhancement of the LIPR due to the short-time dynamics of the 
system. Since we are considering the surface-of-section wavefunctions, the classical dynamics that is important is the 
classical Poincare map Mci. However, since the effects we are seeking are quantum mechanical, we need to start from 
the corresponding quantum Poincare map. For this we use the operator M{k) that enters into the boundary integral 
method: 

{q\M{km= j>dq'K,^{q,q';k){q'\4>) , (24) 



where the kernel is vilA 



KUq.<i; k) ^ ~Y"['\k \v{q) - r(g')l) Hq) ' Hq) - r(<z')], (25) 
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r{q) is a point on the billiard wall at arc length q, and n{q) is the inward-pointing unit normal vector at q. Does M{k) 
as defined generate the correct quantum dynamics? In a sense that is a philosophical question since the Poincare 
map has no obvious exact quantum counterpart. However M{k) satisfies enough properties expected of a quantum 
Poincare map that we shall use it: 

1. To within a semiclassical approximation it agrees with the semiclassical Poincare map (which is easier to 
define 

2. It is approximately unitary within a subspace whose dimension is given by the area of the Poincare section in 
units of Planck's constant, and exponentially small beyond that dimension. 

3. If a boundary function |n) corresponds to a wavefunction that satisfies the boundary conditions at wavenumber 
k = km then \n) is mapped onto itself under the action of M{k): 

M{kn) \n) = \n) . (26) 



Indeed, this is the boundary integral method criterion for an eigenstate. 

The following derivation should be considered suggestive rather than rigorous. For example, M{k) is approximately 
unitary but it has actual eigenvalues both inside and outside the unit circle. This leads to difficulties analytically 
continuing it to the region of interest. 

Given properties 2 and 3 above, it is clear that one crude way to find surface eigenstates of the system would be to 
apply M{k) repeatedly to an arbitrary initial state and average the results: 



lim V[Af(fc„)-eP|0) (27) 
G{kn) !</.) , (28) 



where 



G{k) = lim -^TT (29) 

is similar to a Green's function. Since the component of |(/)) that projects onto \n) is unchanged by the action of 
M{kn), whereas the other components are multiplied by a complex phase, the only component that is not averaged 
out by the above procedure is the projection on |n). 

Whenever k passes through an eigenvalue fc„, one of the eigenvalues e*""''^-' of M{k) passes along the unit circle 
through 1 and G{k) has a singularity. The singular part of G can be written schematically as 

^ Tir.nfi,\ \^ \n)5{k~kn) {n\ 
-ReG{k)=}_^ |rfa„(fc)/dfc| ■ 

The denominator gives the velocity of rotation of the nth eigenvalue as it passes through 1. It can be shown from 
semiclassical arguments that the average of these velocities is given by the average length of the trajectory segments 
corresponding to one Poincare map: 

The second equality follows from a classical theorem well known in acoustics jTsj, where A is the billiard's area 
and L the length of its perimeter. For a chaotic system, Eq. (|3^) holds not only in an average sense but also 
holds approximately for individual eigenphases. Therefore we take {£) outside of the summation, yielding the useful 
expression 

^RcG(fc)«^|n)(5(fc~/c„)(n|. (32) 

n 

Note that the sum is over states of the full system, each with its distinct eigenvalue kn- 
it is a simple warm-up exercise to compute the denominator of the LIPR: 
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^\{<l:P\n)\^ ^^{q.PW) {n\q,p) (33) 

n n 

= {l,P\n) I dkS{k - fc„) {n\q,p) (34) 

= j J dkRe{q,p\G{k)\q,p) (35) 

= J^^J dk[{q,p\q,p) + {q,p\M{k)\q,p) + {q,p\M\k)\q,p) + ■ ■ ■] . (36) 

Now note that M{k) involves phases that vary rapidly with k; therefore only the first term in the brackets survives 
the integral over k. We are left with 

Y^\{q,p\n)f^^{q,p\q,p), (37) 

n 

where Afc is the range of k over which the averaging is done. 
To compute the fourth moment, we proceed as follows: 

X! I(9'-Pl"-)I^ = X! (^iPI") {n\q,p) {q,p\n) {n\q,p) (38) 

11 n 

. /<,M,.-.„)(„|,.rt(,.p|™) /<,r.(t'-.„,(™|,.,).„,„. (39) 

Because of the 5{k — fc„) and 5{k' — km) factors, the Kronecker delta 5nm may be replaced by C5{k — k') where C is 
an undetermined constant with dimensions of momentum which is needed to make its product with the delta function 
dimensionless. Continuing, then, we have 

Y,\{qMn)t ='^{j) J dkdk'Re{q,p\G{k)\q,p)Re{q,p\G{k')\q,p)S{k-k') (40) 

dk[Re{q,p\G{k)\q,p)f (41) 
= c(j) J dkJ2Mq,P\M'{k)\q,p)Re{q,p\M^'{k)\q,p) (42) 

^^(l) J d^'T.[^''(^'P\^^'(''^\l'P)]' (43) 

«f (^)'/ dkJ2\{q,p\A'P{k)\q,p)\' . (44) 

where in going to Eq. ( ^3| ) we have used the fact that for j ^ j', the two phases in the integrand are approximately 
uncorrelated, and the integral averages to zero. Eq. ( ^4[ ) (in obtaining which we have used the fact that for j the 
phase of {q,p\M^ (k)\q,p) is a rapidly changing function of k, and so may be regarded as random in the integration; 
the j = term is irrelevant in the infinite sum) may be interpreted as expressing the proportionality between the 
LIPR and the mean quantum return probability, averaged over all times. As a result of the reloading effect Q the 
mean long-time return probability is proportional to the sum of the short-time recurrences, which may be cut off at 
the mixing time, jmix = Tn-,i^/TB: 

EK9'PI")I'« y(^)V* £ \{q,p\A'P{kUp)\\ (45) 



G 



The constant of proportionality changes in going from Eq. ( [14| ) to Eq. ([45[), as is reflected in our change of notation 
from C to C". The mean quantum return probability is inversely proportional to the number of Planck-sized cells in 
phase space, or inversely proportional to the Heisenberg time, and thus C" is inversely proportional to the Heisenberg 
time. 
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But at times shorter than the mixing time, the quantum map M^{k) can be approximated semiclassically in terms 
of the classical Poincare map M^j (which does not depend on k) times a phase that varies rapidly with k: 



{q,p\M\k)\q,p) « [q,p M^j exp(iHj), 



(46) 



where |q,p) represents a classical Gaussian probability distribution centered at (g,p) (notationally distinguished from 
a quantum wavefunction by the round brackets), is the classical Poincare map iterated j times, and is the 



length of the corresponding classical trajectory, {^q^p 
with the j-times iterated distribution: 



(<?,p|Af^l|g,p) = j dq'dp'P^,^^^^^^^^{q\p')P,^p{q',p') 



q,p] is the classical overlap of the original distribution 



(47) 



So finally we have 



p\n) 



2 



(48) 



Comparing Eq. (p^) and Eq. (Eq) we obtain for the LIPR 



LIPR: 



C'N Ef=lj„_ 










9,^ 





(49) 



(50) 



where the projections in both the numerator and in the denominator are classical. The constant C" can now be 
fixed by requiring that the LIPR, when evaluated far from any short periodic orbits (where scarring plays no role), 
should give the RMT value which we call LIPR(RMT). [As discussed above, LIPR(RMT) = 2 in a system with no 
time-reversal symmetry; in the presence of time- reversal symmetry, LIPR(RMT) reaches the value of 3 on certain 
symmetry lines, as we have seen in Section [VP . The LIPR contribution obtained from short-time dynamics must 
always be multiplied by the appropriate factor given by symmetry considerations; LIPR(RMT) is inherently an effect 
arising from Heisenberg-time behavior and so must be considered separately from short-time contributions to the 
LIPR.] 

Away from short periodic orbits there will be no recurrences, so only the j = term in the above sum contributes. 
Thus, C" = (2Afc/7V) LIPR(RMT) and 



LIPR(g,p) = LIPR(RMT) 





^Cl 









(51) 



The J 7^ terms of the sum give an enhancement to the LIPR whenever a point (g, p) is mapped near itself by an 
iterate of the classical Poincare map; i.e., near periodic orbits. In the next Section we will describe how to evaluate 
the LIPR in the neighborhood of a periodic orbit. 



2. Periodic orbits, linear theory 



The dynamical description of the LIPR developed in the previous Section provides a method of computing the 
LIPR classically, anywhere in phase space, by computing overlaps of classical Gaussian probability distributions with 
themselves iterated under the classical Poincare map Mci- We found that the LIPR is enhanced when the short-time 
overlap is large. This of course is most dramatically the case in the neighborhood of periodic orbits. In this Section 
we show how to compute the LIPR in the neighborhood of a periodic orbit, using only the properties of that one 
orbit. We also generalize to the case of stable and unstable manifolds with arbitrary orientation. 
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Fig. ^ confirms the omnipresence of scarring in the eigenstates of the stadium billia rd, as we no w argue (the 



identification of peaks in Fig. |2| with specific short periodic orbits will be made below in Sec. [VG, IV G 2 in particular). 
In the general theory of scarring, a wavepacket launched on a periodic orbit has a local density of states with a linear 
or short-time envelope (coming from dynamics linearized around the orbit) that consists of an oscillating function of 
energy; isolated homoclinic recurrences, if present, will give rise to further fluctuation multiplying this envelope [Q. 
The concentration of the local density of states in some energy regions and avoidance of others leads to an enhanced 
LIPR in Eq. (|^). Thus the peaks in the LIPR plot in Fig. || at the positions of short periodic orbits must be attributed 
to scarring, even though for the most part the scars are too weak to be visible in coordinate-space plots of individual 
wavef unctions, such as in Fig. |. The scars are there in the sense of non-x^ variation of wavefunction intensity on 
these periodic orbits, as one scans through the eigenstates of the system. Since there are several periodic orbits with 
appreciable scarring, any given eigenvalue is likely to lie near a maximum in the linear envelope of at least some of 
them, and the corresponding eigenstate will more likely than not be scarred on these orbits, and antiscarred on the 
others. Examination of the spectra corresponding to different short periodic orbits reveals that, first, as one varies 
the periodic orbit the position of the smooth envelopes in the local density of states varies so as to cause any given 
energy to lie under the peaks of the envelopes for several different orbits. Secondly, we find that most of the envelopes 
in the local density of states are filled in an egalitarian way (for a discussion of egalitarian versus totalitarian filling 
of the spectral envelopes, see Kaplan and Heller [Q), which increases the likelihood of at least weak scarring along 
that orbit for each eigenstate with energy near the maxima of the envelope. For typical scarred wavefunctions, the 
scarring will not be visible in coordinate space, even though the phase space intensity is greatly enhanced at the 
periodic points, as as one can see in the Husimi plots. Occasionally, the scarring may be strong enough to be also 
visible in a coordinate-space plot; an illustration of strong multiple scarring can be found in the state at A: = 105.608, 
which is simultaneously scarred along the bowtie and diamond orbits; cf. Fig. |l](c). 

In summary, we find a distribution of wavefunction intensities on a periodic orbit which differs from the ^ dis- 
tribution that would be expected on the basis of random-matrix theory, as illustrated in Fig. ^ In Fig. ^(a) we 
show a histogram of wavefunction intensities measured using a Gaussian centered on the bowtie orbit for all states 
lying between fc = 100 and fc = 150. Because the bowtie orbit lies on a symmetry line, the random-matrix-theory 
prediction is a distribution in one degree of freedom (where we have set 6% of the intensities to zero, corresponding 
to the bouncing-ball states), which fails to account for the large number of high- intensity and low-intensity peaks in 
Fig. |^(a). In the presence of scarring, the distribution must be modulated by the linear spectral envelope, which, 
for a Lyapunov exponent of 2.29 for the bowtie orbit, varies from 0.4 to 2.3 periodically in momentum. The period in 
energy is given by 2'k/T where T = 2.6Tb is the period of the bowtie orbit. As a result of this modulation, there are 
many more very small and very large intensities in Fig. ^(a) than would be expected from the naive Porter-Thomas 
(RMT) picture. As plotted on a logarithmic scale, the effect is more pronounced for the larger intensities, but even 
at low intensities we expect an enhancement by about 20%, which is what is indicated by the data between (/)/500 
and (/)/10, where the statistics are good; below (/)/500 the statistics are too poor to distinguish rehably between 
the scar theory and random matrix theory predictions. The amount of strong scarring observed is also quantitatively 
consistent with the predicted scarring corrections to random-matrix theory: for instance, the number of wavefunction 
intensities greater than 10 times the mean is found to be 10 numerically, while the scar-corrected prediction is 13 ± 3 
and the uncorrected random-matrix prediction is only 3±3. At a generic point in phase space, not lying on a scar or on 
a symmetry line, we expect a distribution in two degrees of freedom, and Fig. ^(b) shows the numerical distribution 
of wavefunction intensities at a generic point in phase space, in satisfactory agreement with the random-matrix-theory 
prediction for intensities greater than ^ (/)/10 (without scarring corrections, since there is no scarring at a generic 
point in phase space, but with 6% of the intensities set to zero, namely those corresponding to bouncing-ball states). 
The bouncing-ball states appear of course in the numerical data at intensities I < (/)/10; the total area between the 
data and the curve in this intensity region indeed turns out to be about 6%). Note the smaller number of very small 
and very large intensities as compared to the scarring case. Fig. ^(a). The 2-degrees-of- freedom distribution also falls 
off much faster at large intensities than the 1-degree-of-freedom prediction appropriate for a symmetry line. 

We now turn to the calculation of the heights of the scarring peaks in Fig. || based on a linear expansion around 
periodic orbits. Consider the surfacc-of-scction map in BirkhofF coordinates; we will be interested in the short-time 
linearized dynamics in the vicinity of a periodic orbit The stable and unstable invariant manifolds, defined as 
the locus of points that approach the periodic orbit under infinite iteration of the mapping forwards or backwards 
in time, respectively, will in general have an arbitrary orientation. In the special case that the invariant manifolds 
are aligned with the coordinate axes centered on the periodic orbit of interest — the simplest possibility — we have 
linearized equations of motion q{t) = q(0)e^* and p{t) = p(0)e~^*, where q and p are measured relative to the periodic 
point. An initial Gaussian corresponding to a classical distribution of probability, Po{q,p) = Ce~'^ I" ^ (where 
C = (2/7r)^/^ is a normalization constant) will map after n iterations to 

(g, p)^C exp [-e-^^^q^a^ - e^^^a V] • (52) 
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The overlap with the initial distribution will then be 



J dq dp Pj{q,p)Poiq,p) = ^^J^^^j ' (^3) 

assuming that we are far into the semiclassical regime, i.e. that k is large enough so that the Gaussian wavepacket 
is well contained in the phase space region where the linearized equations of motion are valid. We note as an aside 
that there is no natural way to choose the normalization constant C based on purely classical considerations, since 
classically the return probability of a distribution is not a meaningful concept. Instead, we are implicitly thinking 
of Eq. ( p3[ ) as representing the square of a semi-classical return amplitude of a Gaussian wavepacket, which is a 
well-defined physical quantity, and is given by the square root of Eq. (|5^), times some irrelevant phase associated 
with the action (in units of h) of the periodic orbit. This fixes the normalization constant C above. 

Due to the reloading effect mentioned in the preceding Section (see also [Q) in which returning homoclinic orbits 
contribute in phase, the mean return probability for long times at this periodic orbit, and hence the LIPR, which equals 
a constant times the long-time return probability, will be equal to LIPR(RMT) times the short-time enhancement 
factor 

oo ^ 

Strictly speaking, the sum in Eq. ^ should extend only over times short compared to the mixing time, i.e. \j\ < 
X^^ log(fcL); however this cutoff is irrelevant in the semiclassical limit kL 3> 1 in which we are working. 

It is desirable to extend this result to the case when the invariant manifolds are not perpendicular. This is equivalent 
to considering perpendicular invariant manifolds but with a non-circular Gaussian (as one can go from one case to 
the other via a canonical transformation). Let x — {q,p) and A be the quadratic form defining the initial Gaussian 
centered on a periodic orbit, 

Po(x) = Ce-^"^^ , (55) 
with det A = 1. We define J to be the Jacobian of A/ci at the periodic orbit; then 

Pj{x) = Cexp [-x^{J-^fAJ-'x] (56) 

and the overlap is 



dq dp P,{q,p)Po{q,p) = 2^ det[A + ( J-^)^A J-^] ^^^^ 

(again, we see that using our normalization the overlap is unity at time j = 0). Thus, a better prediction of the 
LIPR peak heights at the periodic orbits — one which takes into account the local stable and unstable directions at 
the periodic orbit, and not only the stability exponent — is LIPR(RMT) multiplied by 



= ^2^Jpet[A+iJ-^)TAJ-^- 

Note that S2 reduces to 5*1 when the invariant manifolds are perpendicular, and the Gaussian test state is oriented 
along these manifolds, as an easy calculation shows. Of course, using a canonical transformation we may easily see 
that for any orientation of the manifolds there is an infinite one-parameter family of 'optimal' Gaussian test states, 
all of which display maximum possible scarring in accordance with Eq. ( |5^ ) . Both Si and 5*2 approach unity for long 
or highly unstable orbits, when the instability exponent A becomes large. 

We have calculated Si and ^2 for theperiodic orbits appearing in Fig. 0, a nd c ompare these theoretical predictions 



with the actual peak heights in Table III. All quantities appearing in Table |II| are either predicted or actual values 
of LIPR(g,p)/2 at a periodic point {q,p)- In other words, we have divided out by the RMT prediction LIPR = 2, 
which is relevant for all of space space away from the symmetry lines. However, on the symmetry lines the RMT 



expectation is 3 rather than 2 (see Section IV D| ), so the theoretical values Si and 5*2 have been multiplied by 3/2 for 



those periodic points which do lie on symmetry lines, in order to compare with the actual LIPR/2 data in the fifth 



and sixth columns of Table III. In fact, most of the points listed in Table [II do lie on one or another symmetry line; 



the first column of the table indicates which symmetry class each point lies in (see definition of the symbols in the 
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Orbit 


Symmetry 


Lyapunov exponent 
A 


Si 


5*2 


actual 
k = 100-150 


height 
k = 200-225 


classical 


HB 


a 


1.76 


2.72 


2.15 


1.86 


2.01 


2.78 


Vi 


b 


2.02 


2.40 


2.36 


1.91 


1.93 


2.36 


Va 


c 




2.40 


1.95 


1.68 


1.87 


1.97 


Di 


a 


2.06 


2.37 


1.97 


1.89 


1.80 


2.03 


Da 


b 




2.37 


2.27 


2.15 


2.06 


2.27 


Bi 


d 


2.29 


2.18 


2.06 


1.86 


1.87 


2.02 


Ba 


d 




2.18 


2.06 


1.86 


1.87 


2.02 


Zi 


e 


2.38 


2.10 


2.04 


1.69 


1.88 


2.04 


Za 


e 




2.10 


2.04 


1.69 


1.88 


2.04 


Z3 


c 




2.10 


1.77 


1.53 


1.55 


1.85 


Xi 


df 


2.58 


2.00 


1.94 


1.91 


2.03 


2.07 


Xa 


df 




2.00 


1.94 


1.91 


2.03 


2.07 


Ai 


b 


2.60 


1.98 


1.86 


1.53 


1.55 


1.85 


Aa 


d 




1.98 


1.86 


1.59 


1.51 


1.90 


A3 


d 




1.98 


1.86 


1.59 


1.51 


1.90 


Wi 


b 


2.68 


1.94 


1.94 


1.71 


1.86 


1.95 


Wa 






1.29 


1.27 


1.68 


1.78 


1.29 


W3 


c 




1.94 


1.70 


1.66 


1.81 


1.84 


W4 






1.29 


1.27 


1.68 


1.78 


1.29 


Hi 


a 


3.26 


1.74 


1.73 


2.01 


2.37 


2.00 


Ha 


f 




1.74 


1.74 


1.78 


2.12 


1.99 


H3 


f 




1.74 


1.73 


1.78 


2.12 


2.00 


Ti 


b 


3.83 


1.64 


1.64 


1.74 


1.89 


1.77 


Ta 


f 




1.64 


1.56 


1.15 


1.26 


1.61 


T3 


f 




1.64 


1.56 


1.15 


1.26 


1.61 


Fi 


a 


4.90 


1.55 


1.53 


1.58 


1.63 


1.95 


Fa 






1.03 


1.04 


1.14 


1.18 


1.02 


F3 


d 




1.55 


1.53 


1.50 


1.62 


1.57 


F4 


d 




1.55 


1.55 


1.50 


1.62 


1.57 


F5 






1.03 


1.03 


1.14 


1.18 


1.02 



TABLE HI. Scar peak heights (values of LIPR/2) on short periodic orbits, as predicted by the linearized semiclassical 
measures Si, S2 and by the exact brute-force classical simulation, defined in the text, and to be compared with actual 
quantum-mechanically computed heights. Periodic orbit labels are defined in Table The symmetry family to which each 
point on each periodic orbit belongs is labeled as follows: a, orbit hits center of the circular endcap; b, hits center of the straight 
segment; c, incident normal to the boundary; d, possesses a vertical segment in the endcap region; e, passes through the center 
of the stadium; and /, possesses a horizontal segment. The Lyapunov exponent is listed only once for each orbit. Again, the 
predictions Si and S2 are based on the tangent map in the desymmetrized stadium. The predictions in columns 3, 4 and 7 
include the quantum symmetry correction factor of 3/2 for periodic points with special symmetry, as indicated in column 1. 
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caption). As expected, 5*1 fails to be an accurate measure of the LIPR peak heights because the test states in general 
are not aligned optimally with the invariant manifolds. With the inclusion of symmetry corrections, the agreement 
of 5*2 is roughly correct, as can be seen in a scatter plot of the predicted versus calculated LIPR/2 peak intensities, 
Fig.]|; as shown in Fig. ^(c), the scatter is markedly worse when the symmetry corrections are omitted. 

In order to appreciate the significance of these results, we should d iscu ss the factors contributing to the uncertainty 



in the measured LIPR/2 values listed in columns 6 and 7 of Table |III| . There are two issues: first, there may be a 



non-integral number of oscillations of the local density of states in the energy window considered. The size of this 
effect scales inversely with the number of scar oscillations in the energy window, and is only a few percent for the 



data presented in Table [II. A second, and bigger, source of error arises from the finite number of eigenvalues in the 
window. For random fluctuations, the uncertainty coming from this effect scales as the square root of the mean level 
spacing divided by the size of the window. This is why the windows from k = 100 to 150 and from k — 200 to 225 were 
taken to be as large as possible. A monte carlo simulation was done of random wavefunction intensities constrained 
by the "linear" spectral envelope: it suggests that the statistical uncertainty in the LIPRs should range from 5% to 
8% (depending on the stability matrix of the orbit in question) for periodic points on symmetry lines, and should be 
slightly smaller for those points that are not on any symmetry line. This level of fluctuation is roughly consistent with 
the dispersion in the scatter plots of Fig. and does not alter the conclusion that 5*2 together with symmetry factors 
is a good prediction of the peak heights. Thus, the symmetry lines emerge as essential to a quantitative understanding 
of the calculated LIPRs due to scarring. A final consideration supporting the validity of our symmetry analysis is the 
pattern of peak heights for the F orbit, on which periodic points Fi, F3 and F4 have a symmetry correction while F2 
and F5 do not; thus, the predicted peak heights for the points having symmetry corrections are about fifty percent 
higher than for the others, and this predicted pattern is reproduced in the calculated heights. 

The brute-force classical calculation at times large compared to Tg goes beyond Si and S2 to include the effects 
of nonlinear homoclinic recurrences, which add to the predicted peak height. For each periodic point, the returning 
probability for an initial Gaussian distribution centered at that point was integrated up to a cutoff time of 5Tb, which 
was chosen to be large compared to Tb but smaller than the mixing time Tmix ~ 7Tb, after which the returning 
probability approaches a constant per unit time independent of position in phase space. In the semiclassical limit, 
where strong, identifiable recurrences at times between the initial decay of the wavepacket and the Heisenberg time 
can be neglected, the classical simulation reproduces 82- The classical simulation reproduces, as it must, one feature 
of the quantum data: the peak heights it predicts agree for points related by symmetry, namely Bi and B2, Zi and 
Z2, Xi and X2, Ai and A2, W2 and W4, H2 and H3, T2 and T3 and F3 and F4. In predicting the quantum heights, 
the brute-force classical calculation does no better than ^2 and makes predictions that are systematically high (after 
including symmetry effects). In order to investigate the energy dependence of the classically predicted peak heights, 
the classical simulation was run both at fc = 100 and at = 200, the difference being in the size of the Gaussian 
as a fraction of the phase space area. For almost all of the periodic orbits the results at the two energies agreed to 
within a few percent, which is of the same size as the discrepancy between the numerical quantum peak heights in 
those two energy ranges. This tells us that isolated, identifiable non-linear recurrences beyond the initial decay of the 
wavepacket away from the periodic orbit are not that important for computing the mean long-time return probability. 
Any such nonlinear effects would be strongly dependent on the size of the initial wavepacket, i.e. on the wavevector 
k. 

The classical calculation for the average return probability can be performed anywhere in phase space, not only on 
the short periodic orbits considered in this Section; see Sec. [V G 3 below. There we shall see that we get a picture 
similar to Fig. 0, but with more diffuse peaks and the onset of mixing for T > Tj^ix- 



F. Dynamics plus symmetry 

In the previous two Sections, we have seen that discrete symmetries have two important and distinct effects on 
wavefunction localization properties. First, in the presence of time-reversal symmetry, the horizontal line of zero 
parallel momentum in the surface of section phase space, as well as vertical lines associated with parity symmetry, 
if any, both have a local IPR of 3, in contrast with the background LIPR of 2 that is seen away from the symmetry 
lines. The values LIPR — 3 and LIPR — 2 correspond to the random matrix theory predictions for real and complex 
Gaussian random wavefunctions, respectively. Furthermore, short-time iterates of these symmetry lines result in 
additional (curved) symmetry lines where the LIPR takes values intermediate between 2 and 3. When a periodic 
point happens to be located on a symmetry line (which in a highly symmetric system such as the stadium happens 
quite often for the short orbits), the LIPR enhancement associated with scar effects must be combined with the 
enhancement due to symmetry. In the semiclassical limit where our statistical analysis is relevant, the two effects are 
independent of one another and simply multiply, the scar effect being associated with short-time dynamics around a 
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periodic orbit, while the symmetry effect appears near the much larger Heisenberg time, where individual eigenlevels 
are resolved. 

Secondly, in the presence of spatial symmetries such as parity, symmetry must be taken into consideration when 
computing periodic orbit properties (period, classical action and monodromy matrix) for the purpose of quantifying 
the scar effect. This holds true even if the system is not time-reversal invariant, so that no symmetry lines are present 
in the LIPR plot. In the previous Section, we have addressed this second issue simply by considering orbits in a 
desymmetrized version of the billiard. That approach, however, is neither rigorous nor completely general, and in 
particular is problematic for orbits like the horizontal bounce, which does not exist in the interior of the quarter- 
stadium. In this Section we present an alternative approach to including spatial symmetry in computing scar effects 
on localization. 

In the stadium there are two parity operations, reflection in x and in y, which we represent by the operators 
and i?y, respectively, and which separate the stadium wavefunctions into four distinct symmetry classes. The sums 
on the left-hand side in Eqs. ( p7| ) and ( ^8| ) can be divided into sums over each symmetry class separately; thus, for 
instance, the left-hand side of Eq. (|3^) becomes 

(59) 
(60) 



n 71-1-4- n_-|- n-|-_ n 

= ^ [\{qMP++\n)\' + \{qMP-+\n)\^ + \{qMP+-\n)? + \{qMP-\n)Y 



where P±± is the projector onto the subspace with even or odd symmetry under Rx and Ry^ respectively: P±± = 
(1 ± i?a;)(l ± Ry)/^- A similar development of the left-hand side of Eq. ( [48|) for the fourth moment is possible. The 
projection operators may now be considered to act to their left, on the state ho, g). This means simply that |p, (?) 
is to be replaced by g), etc., on the right-hand sides of Eqs. (|3^) and (^sj). If we consider the LIPR for the 

even-even states only we obtain 



LIPR++(g,p) 



N++Y.n\{iMP++\n)\' 



Y.n\{<lMP++W) 



LIPR(RMT)- 



LIPR(RMT) ^ 



LIPR(RMT) ^ 



(g,p|P++|9,p)P 



{q,p\P++\q,p) |2 



4 J2ii,i2=0.1 (^'P 



q,p 



1/2 



{<i,p\P++\q,p) 



(61) 

(62) 
(63) 

(64) 



The constant of proportionality here, C'_^^ = (2Afc/7V+_|_)LIPR(RMT), is four times larger than C in Eq. (p5| ) because 
A''_|_+ = N/A] this reflects the fact that the Heisenberg time is four times smaller for the desymmetrized stadium (C;^+ 
being inversely proportional to Th++)- 

Thus we see that in addition to the usual recurrences due to periodic orbits {A'P {q,p) = {q,p), ii, 12 = 0) there are, 
for certain initial values {q,p), recurrences when the orbit closes up to a reflection in x or y or both. For example, the 
V orbit starts at the center of the straight segment, hits the semicircular segment at j = 1, and at j = 2 returns to 
where it started up to reflection in x. Another example is the Z orbit, which passes through the center of the stadium; 
at J = 3 it returns to where it started up to a reflection in both x and y. 

Similar recurrences up to a symmetry operation would be present in the LIPR summed over all symmetry classes. 
One last important point to note is that for an orbit located right on a symmetry line (such as the horizontal bounce 
orbit in the stadium), the corresponding symmetry {Ry in this case) is not included in computing the scar effect. 
The fact that the wavepacket is located on a symmetry line is of course relevant in setting LIPR(RMT)= 3 instead 
of 2. The other symmetry, namely R^, still needs to be included: for the horizontal bounce orbit its effect is to make 
short-time recurrences in Eq. ^ happen for all values of j instead of for even j only. 
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G. Derivation as return probability 



Finally, we explain the interpretation of the LIPR as the infinite-time average return probability, as in Eq. ( |65| ) 
below (where it is assumed that the energy window in the sums over n and n' is large enough to cover the energy 
scales present in the test state \a)). In order to clarify what is happening in the infinite-time limit it is useful to 
look at finite times, and then let time tend to infinity. Thus we define the average return probability at time T for a 
wavepacket to be |lq| 



1 

2T 



T 

dt 

-T 

T 



J ^dt\{^a\n)\'\{^a\n')\\ 



In the limit as T tends to infinity. 



sin((£:„-£:„Or) 

{En — En')T 



5{E.n ~ E„ 



and (in the case of a non-degenerate spectrum) we recover Eq. (|^) up to an overall proportionality constant. The 
arbitrary state |^'a) living in the interior of the stadium may be replaced by a Gaussian \q,p) living on the boundary 
if at the same time the eigenstatcs |5'„) are replaced by their normal derivatives on the boundary, \n). The natural 
timescale for T is the time Tb for bouncing between the straight segments of the billiard. With k ~ 100 and m = 1/2 we 
have V = p/m — Tik/ni = 200 and so Tb — 2/w — 0.01. For T Tb we should have Paa{T) ^ l/T; for intermediate 
times we expect classical behavior as long as T < Tmix, where the mixing time Tmix is the time it takes for the 
classical dynamics to spread a given phase-space element throughout the entire phase space (in order of magnitude it 
is given by Tb times the logarithm of the number of Planck-sized cells in the classical surface-of-section phase space, or 
~ Tb log(fcL) ). The classical behavior will be computed directly below, where based on the time needed for a simulated 
classical ensemble to begin to spread smoothly throughout phase space we find Tmix ~ 7Tb at /c = 100; that it is so large 
may be attributed to intermittency due to bouncing-ball orbits. Finally for T ^ Tb-, that is beyond the Heisenberg 
time where individual eigenstates are resolved (here, Th = 1/AE = l/{2kAk)^ 1/(2 • 100 • 0.05) = 0.1 = lOTg at 
k ~ 100), we should find that Paa{T) tends to the LIPR plot as shown in Fig. ||. 

1. Finite-time return probability 

The results for the time-dependent average return probability P{q,p),(q,p) (^) are given in Fig. ^, for the energy window 
at k — 100. We are unable to study numerically the limit T <C Tb because of the finite number of eigenfunctions 
included in the sum in Eq. (|65|). At T = Tb, however, we see large amplitudes coming from the bouncing-ball and 
whispering-gallery regions of phase space, both of which are short-time classical effects. At T = 2Tb the horizontal- 
bounce orbit appears. By T — 3Tb several discrete peaks appear, and more of them are present at T = 4Tb. These will 
be shown below to correspond to short periodic orbits of the classical dynamics in the desymmetrized quarter-stadium 
billiard. As T becomes as large as SiTb and then IOTb the finite-time average return probability indeed approaches the 
infinite-time LIPR. The non-uniform structure of the infinite-time return probability is a purely quantum-mechanical 
effect, because according to classical mechanics the average return probability must become uniform for T 3> Tmix- 

2. Identification of peaks 

Let us work with the data at T = 4Tb, where the strong peaks are all present and readily distinguishable from 
the background [see Fig. ||(d)]. We identify these peaks with the short periodic orbits in Table IV and Fig. 0. To 



avoid having to deal explicitly with symmetry effects, we use the properties of the reduced orbits in the fundamental 
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Label 


Description 


q 


p/k 


Length 
(dcsymmetrized) 


Period/Ts 
(dcsymmetrized) 


HB 


horizontal bounce 


2.5708 





1 


2.00 


V 


V-shaped 





0.70711 


2 


2.41 


D 


diamond 


2.5708 


0.44721 


2 


2.24 


B 


bowtie 


1.5236 


0.5 


2 


2.60 


z 


Z-shaped 


0.5 


-0.44721 


3 


3.24 


X 


box 


1.7854 


0.70711 


2 


2.41 


A 


accordion 





0.5491 


3 


3.30 


W 


W-shaped 





0.31623 


4 


4.16 


H 


hexagon 


2.5708 


-0.86603 


3 


2.50 


T 


triangle 





0.83573 


3 


4.30 


F 


fish 


2.5708 


0.6478 


5 


5.40 



TABLE IV. Principal short periodic orbits appearing in Figs. ^ and |9|. The initial coordinate q is the distance around the 
perimeter of the stadium measured from the center of the upper straight segment; the initial momentum p is positive in the 
clockwise sense. The lengths and periods are given for the desymmetrized stadium, which is relevant to our calculations. 



domain — one quarter of the stadium. (This approach has been justified rigorously in Sec. IV F.) Thus the periods 
of the most important orbits, as tabulated in Table IV, range from 2Tb to 5Tb although their periods in the full 
stadium may be longer. Although the triangle is a short periodic orbit, it does not appear until 5Tb because its 
desymmetrized version is just as long as the orbit in the full stadium. As T increases beyond Th ~ lOTs, many of 
the short periodic orbits become incorporated into the streaks that dominate the IPR plots. 



3. Classical finite-time return probability 



The persistence to long times of short-time classical information is a remarkable property of the quantum mechanics. 
Classically, we would expect the recurrences to be washed out at times large compared to the mixing time. Indeed, we 
have computed a classical analogue to the time-dependent average return probability, namely the overlap of an initial 
Gaussian in classical (Birkhoff) phase space with its iterates under the classical mapping. In order to get results in 
continuous time, the overlap of the iterate at (real) time t with the initial Gaussian was stored as a delta function at 
time t and at the end of the calculation the integral over t from — T to +T was taken, as in the quantum calculation. 
As shown in Fig. |l^, by T = lOTg all of the periodic orbits seen in the quantum calculation are present, although 
the peaks are more diffuse classically, but by 20Tb the classical picture is beginning to be uniformly distributed 
throughout phase space, losing memory of the short-time behavior. Quantum-mechanically, of course, the peaks do 
not become washed out at long times and persist in the infinite-time limit, Fig. 



V. CONCLUSION 



The advantage of the approach implemented in Section |^ is that all of the scars can be observed simultaneously 
in the plot of inverse participation ratio as a function of phase space location. Although scars are associated with 
the classical structure of unstable periodic orbits, the persistence to long times of an excess in the average return 
probability is a purely quantum-mechanical phenomenon. The brightnesses of the peaks in the LIPR plot can be 
predicted semiclassically with some success (to accuracy of about 10 — 15%), the size of the deviations being roughly 
consistent with expected statistical fluctuations. More precise comparison with the scar theory predictions would 
require producing a greatly increased set of eigenstates: this could be done by going to higher energies or (with far 
less computational expense) by using an ensemble of chaotic systems with a given unstable orbit k ept fixed 0. The 



were no more 



predictions of brightnesses based on exact classical phase-space evolution at the end of Section |V 
successful than those based on linearized behavior around the periodic orbit. 

Most of the extra localization of individual eigenfunctions in Husimi phase space, over and above Gaussian random 
fiuctuations, may be satisfactorily accounted for by quantum symmetry effects. Scarring is also seen to be present 
based on an analysis of the spectra of wavepackets located on short periodic orbits. Symmetry effects on the overall 
mean IPR scale as 0{Vh) (i.e., 0{y^)), while the effect of scarring on the overall mean IPR scales as 0{h) (i.e., 
0(l/fc)). The reason for this, for the symmetry effects, is that the widths in position and momentum of the test-state 
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Gaussian scale as u ~ \ffi^ setting the width of the symmetry Hnes (Eqs. (]20|), (|2l|)), while the length of the symmetry- 
enhanced lines is 7i-independent. Thus, the total fraction of phase space covered by the symmetry lines scales as \ffi. 
The scar effect, on the other hand, is significant for those phase-space Gaussian wavepackets that have large intensity 
on the periodic orbit and the fraction of such Gaussians goes as the area of the Gaussian; i.e., as \/h, x = U. 
We note that although the size of both these effects on the overall wavefunction IPR is expected to go to zero in 
the semiclassical limit (because periodic orbits and symmetry lines both affect a smaller and small fraction of phase 
space surrounding them, in the ?i — > limit), the local IPR (LIPR) either on a symmetry line or on a periodic orbit is 
?i- independent. Thus, as measured using the local IPR, deviations from naive RMT predictions persist to arbitrarily 
high energies and do not decay in the semiclassical limit. As for other possible kinds of eigenfunction localization, 
not associated with the short periodic orbits of the system or with symmetry lines, their explanation remains an open 
question. 
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FIG. 1. Representative eigenfunctions of the stadium billiard in coordinate space (left) with x-axis running from to 4 and 
y-axis running from —1 to 1, and in Husimi phase space (right) with distance q along the perimeter (horizontal axis) running 
clockwise from (corresponding to the center of the upper straight segment) to 4 + 27r and tangential momentum (vertical axis) 
running from to where k is the wavenumber of the eigenfunction. In the coordinate space plots, the graylevel represents 
\ip\'^. In the Husimi phase space plots, white is high intensity and black low intensity, (a) k = 24.680, (b) k = 100.787 and (c) 
k = 105.608, chosen to illustrate scarring along the bowtie and diamond orbits (cf. Table IV). 
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(c) 



FIG. 2. LIPR (see Eq. in phase space as a function of q running horizontally from to 4 + 27r and p/k running vertically 
from —1 to 1. The test-state width a is chosen so that the test Gaussians have circular cross-sections when the figures are 
plotted with a square aspect ratio. Darker shading indicates higher values, (a) fc = 50 to 60, (b) k — 100 to 150, and (c) 
k = 200 to 225. 
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FIG. 3. LIPR in phase space for eigenfunctions in the range from k = 100 to 150, coordinates as in Fig. |2|(b), for (a) 
~ 0.01, position-like test Gaussians, to be compared with ~ 0.05 used in Fig. ^(b), and (b) ~ 0.25, momentmn-like 
test Gaussians. 
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FIG. 4. Principal symmetry-related structures in phase space of the inverse-participation ratio plotted in Fig. ^ and the 
average-return probability plotted in Fig. ^. Schematic depiction of the symmetry lines in Fig. ^, labeled in the order in which 
they are described in the text. 




FIG. 5. LIPR for real test wavepcakets (see Eq. (^2|)) in phase space as a function of q running horizontally from to 4 + 27r 
and p/k running vertically from —1 to 1. The axes and parameters are the same as those used in Fig. ^ 
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FIG. 6. Probability distribution of wavcfunction intensities In ~ \{q,p\n)\ for the 1746 energy levels between k = 100 
and 150. (a) For {q,p) situated on the bowtie orbit. Solid line, numerical distribution of wavefunction intensities; dashed 
line, scar-corrected prediction of random-matrix theory with one degree of freedom; dotted line, random-matrix prediction, 
uncorrected for scarring. For purposes of illustration, the low-intensity part is shown in the figure on the left in a semilog plot 
and the high- intensity tail on the right in a log-log plot, (b) Same as (a) but for {q,p) at a generic point in phase space, and 
the dashed line is now the random-matrix theory prediction for two degrees of freedom (without scar corrections, but adjusted 
for the 6% nearly vanishing intensities due to bouncing-ball states). 
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FIG. 7. Principal dynamics-related structures in phase space of the inverse-participation ratio plotted in Fig. g and the 
average-return probability plotted in Fig. ^ The labeling is according to Table |^ notation of scars on the principal short 
periodic orbits in Fig. ^(e) , in an irreducible one-eighth region of phase space. The other regions of phase space are obtained 
by symmetry. The most prominent scarred orbits are W, Z, D, B, X and H. 
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FIG. 8. Scatter plots of computed peak heights in columns 5 and 6 of Table [II versus the linearized semiclassical predictions 
S2 in column 4 and the brute-force semiclassical prediction in column 7. (a) Eigenfunctions ranging from fc = 100 to 150, 
symmetry corrections included in 52; (b) k — 200 to 225, symmetry corrections again included in 5*2; (c) k = 100 to 150, but 
symmetry corrections not included in 82- Without symmetry corrections the agreement is much worse, (d) k — 100 to 150 
versus brute-force classical prediction without symmetry correction and (e) with symmetry corrections included. 
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FIG. 9. Quantum-mechanical average-return probability PaaiT) as defined in Eq. (|65|), near k = 100; g and p as before, (a) 
T = Tb, (h) T = 2Tb, (c) T = STs, (d) T = 4rs, (e) T = STs and (f) T = lOTs. Note the emergence of periodic orbits at 
T = 2Tb to T = STs. For T = 10Tb the figure is already beginning to resemble the LIPR plotted in Fig. |(b). 
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FIG. 10. Results of a classical simulation of the return probability versus time, q running horizontally from to 4 + 27r and 
p vertically from —1 to 1. Initial Gaussians were chosen with = 0.05, corresponding to a spread in position space of 0.2 out 
of 4 + 27r, or 2% of the perimeter. All of the major periodic orbits are present here as short-time recurrences, (a) T = 5Tb, (b) 
r = IOTb, (c) T = 2QTb. 
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